function [R,F,Ab,LK] = resp2dexample(keff,order,mat)
% this function will produce the response matrix for
% the simple 2-D problem given a keff and legendre order
%keff=1;order=1;mat=1;
global numg numel
it = 0;
%        D       sigR   chi*nu*sig   scattering
dat = [ 1.500   0.030   0.000        0.000  0.000  
        0.420   0.080   0.125        0.020  0.000
        2.000   0.040   0.000        0.000  0.000
        0.300   0.010   0.000        0.040  0.000 ];    
abdat = [0.010
         0.080
         0.000
         0.010];

numm = 2;

if numel == 1
    % simple one coarse mesh element
    Y=10;
    xcm  = [   0     2     4      6  ];
    xfm  = [     1      1      1     ]*Y;
    ycm  = [   0     2     4      6  ];
    yfm  = [     1      1      1     ]*Y;
mt   = [  1 1 1 
          1 2 1 
          1 1 1  ]';  
else
    Y = 10;
    % still simple, but nine coarse meshes, uniform
%     xcm  = [   0     2 ];
%     xfm  = [      2    ]*Y;
%     ycm  = [   0     2 ];
%     yfm  = [      2    ]*Y;
%     mt   = [  mat  ]; 
    % simple one coarse mesh element
    xcm  = [   0     2     4      6  ];
    xfm  = [     1      1      1     ]*Y;
    ycm  = [   0     2     4      6  ];
    yfm  = [     1      1      1     ]*Y;
    mt   = [  1 1 1 
              1 2 1 
              1 1 1  ]';
end

BCL = 0;
BCR = 0;
BCB = 2;
BCT = 0;      

src(:,:,1) = [ 0 
               0 ];            
src(:,:,2) = [ 0
               0 ];           

IBSL = [ 0 
         0 ]';
IBSR = [ 0
         0 ]';
IBSB = [ 1 
         0 ]';
IBST = [ 0 
         0 ]';
%  4 faces * groups * orders 2*4*2
R = zeros(4*numg*(order+1));
Rinc = zeros( numg*(order+1) );
Rfar = Rinc; Rlef=Rinc; Rrig=Rinc;
F = zeros( 1, numg*(order+1) );

% loop over legendre order and group --------------------------------------
for RG = 1:numg     
for RO = 0:order

% -------------------------------------------------------------------------
[AB, AL, AC, AR, AT, S, dx, dy, dv, dc, siga, nusig, sct] = ...
    twoDcoefMGresp( dat, abdat, numg, numm, xcm, xfm, ycm, yfm, mt, src, it, ...
	[BCL BCR BCB BCT], RO, RG );

N = sum(xfm); M = sum(yfm);
for i = 1:N+1
    for j = 1:M+1
        k = i+(j-1)*( N+1 );
        indx(k,1)=i;
        indx(k,2)=j;
    end
end

xint = zeros(N+1,1); 
for i = 2:N+1
    xint(i)=xint(i-1)+dx(i-1);
end
% pad end of AB,AL
AB = [AB' zeros(numg,N+1)]';
AL = [AL' zeros(numg,1)]';
% pad top of AR,AT
AT = [zeros(numg,N+1) AT']';
AR = [zeros(numg,1) AR']';
phi = zeros( (N+1)*(M+1), numg );
scsrc = zeros( (N+1)*(M+1), 1);
%==========================================================================
% set A1 and A2
A1 = spdiags([AB(:,1) AL(:,1) AC(:,1) AR(:,1) AT(:,1)], ...
    [-N-1 -1 0 1 N+1],(N+1)*(M+1),(N+1)*(M+1));
A2 = spdiags([AB(:,2) AL(:,2) AC(:,2) AR(:,2) AT(:,2)], ...
    [-N-1 -1 0 1 N+1],(N+1)*(M+1),(N+1)*(M+1));
% solve the first group flux
phieps = 1e-8;  mxit = 60;  it = 0; phierr = 1;

while (phierr > phieps && it < mxit)
phi0 = phi;
%A = spdiags([AB(:,1) AL(:,1) AC(:,1) AR(:,1) AT(:,1)], ...
%    [-N-1 -1 0 1 N+1],(N+1)*(M+1),(N+1)*(M+1));
phi(:,1) = A1 \ ( S(:,1) +  1/keff * nusig(:,2).*phi(:,2) );
for g = 2:numg
    scsrc   = zeros( (N+1)*(M+1), 1);
    for gg = 1:g-1
        scsrc(:,1) = scsrc(:,1) + sct(:,g,gg).*phi(:,gg);
    end
    %A = spdiags([AB(:,g) AL(:,g) AC(:,g) AR(:,g) AT(:,g)], ...
    %    [-N-1 -1 0 1 N+1],(N+1)*(M+1),(N+1)*(M+1));
    phi(:,g) = A2 \ ( S(:,g) + scsrc(:,1) );
end
phierr = norm( phi0 - phi );
if max(max(phi)) > 100
    disp('somewhere, phi has blown up!')
end
it = it + 1;
end
if (it>=mxit)
    disp('exceeded mxit...')
end
%==========================================================================

% responses of various sorts, reponse due to  RG,  RO incident current
[CRL CRR CRB CRT LL LR LB LT] = cresp(phi,dx,dy,dc,mt,xfm,yfm,numg);
f1 = sum( phi(:,1).*nusig(:,1) );
f2 = sum( phi(:,2).*nusig(:,2) );
F(1, (RO)*numg+RG) = f1+f2;
a1 = sum( phi(:,1).*siga(:,1) );
a2 = sum( phi(:,2).*siga(:,2) );
Ab(1, (RO)*numg+RG) = a1+a2;
% compute leakage
LK(1,(RO)*numg+RG) =  LB; %incident
LK(2,(RO)*numg+RG) =  LT; %far
LK(3,(RO)*numg+RG) =  LR; %to right of incident
LK(4,(RO)*numg+RG) =  LL; %to left of incident

% determine legendre coefficients -----------------------------------------
N = sum(xfm);    % max degree
k = 0:N;         % interval
maxord = order;  % maxorder to which we compute P's
P = zeros(N+1,N+1);
P(1,1:N+1) = 1; 
P(2,1:N+1) = 1 - (2*k)/N;
for j = 3:N+1 %otherwise N-1
    i = j-2;
    P(i+2,1:N+1) = ( (2*i+1)*(N-2*k).*P(i+1,:) - i*(N+i+1)*P(i,:) ) ./ ...
        ((i+1)*(N-i)) ;
end
for i = 0:maxord 
    a = (2*i+1);
    b = factorial(N+i+1)/factorial(N);
    c = factorial(N)/factorial(N-i);
    w(i+1)=a*c/b;
end
x = xint;

for G=1:numg
    % f def'n ASSUMES four-way symmetry
    f(1,:) = CRB(end:-1:1,G)';
    f(2,:) = CRT(:,G)';  
    f(3,:) = CRR(end:-1:1,G)';
    f(4,:) = CRL(:,G)';
 
    % now we want to expand f into fi
    for ord=1:maxord+1
        ff(:,ord) = w(ord) * f(:,:) * P(ord,:)';
    end
    side=1;
    %figure(10)
    %fplt = zeros(1,N+1);
    %for ord=1:order+1
    %    fplt = fplt+ff(side,ord)*P(ord,:);
    %end
    %plot(k,f(side,:),'k',k,fplt,'r.')
    %legend('f','order',0)
    %G:numg:4,  (order+1)*(RO)+RG
    Rinc( G:numg:end, (RO)*numg+RG ) = ff(1,:)';
    Rfar( G:numg:end, (RO)*numg+RG ) = ff(2,:)';
    Rrig( G:numg:end, (RO)*numg+RG ) = ff(3,:)';
    Rlef( G:numg:end, (RO)*numg+RG ) = ff(4,:)';
   % Rinc=abs(Rinc);
end

% end loop over legendre order and group ----------------------------------
end
end
% -------------------------------------------------------------------------
R = [  Rinc  Rfar  Rrig  Rlef
       Rfar  Rinc  Rlef  Rrig
       Rlef  Rrig  Rinc  Rfar
       Rrig  Rlef  Rfar  Rinc ];               
F = [ F F F F];
Ab= [Ab Ab Ab Ab];
